# ===============================================================#
#                     Replication files for:                     #
#.  "Attitudinal and Behavioral Legacies of Wartime Violence:    #
#                      A Meta-Analysis"                          #
#                        Joan Barceló                            #
#               American Political Science Review                #
#               Last update: September 3, 2025                   #
# ===============================================================#

#################################
# Figure 2: Figure on Engagement, Trust and Prosociality Outcomes
################################

## ------------------------ Load Packages ------------------------

library(tidyverse)
library(patchwork)
library(grid)

# ---------- Paths ----------

base_path_meta    <- "~/Datasets/"
base_path_effects <- "~/Results/"

read_if <- function(path) if (file.exists(path)) read.csv(path) else NULL

# Styles
col_full <- "indianred3"
col_biv  <- "steelblue4"
col_qe   <- "forestgreen"
sh_full  <- 15
sh_biv   <- 15
sh_qe    <- 17
lbl <- function(b, p) paste0(round(b, 2), ifelse(p < .001, "***", ifelse(p < .01, "**", ifelse(p < .05, "*", ""))))

## ------------------------ Load datasets ------------------------
# Groups
meta.groups               <- read_if(paste0(base_path_meta, "meta.groups.csv"))
meta.groups.random        <- read_if(paste0(base_path_meta, "meta.groups_random.csv"))
meta.groups_biv           <- read_if(paste0(base_path_meta, "meta.groups_biv.csv"))
groups.effects            <- read_if(paste0(base_path_effects, "7.results_group_participation_models.csv"))

# Participation
meta.participation        <- read_if(paste0(base_path_meta, "meta.political_participation.csv"))
meta.participation.random <- read_if(paste0(base_path_meta, "meta.political_participation_random.csv"))
meta.participation_biv    <- read_if(paste0(base_path_meta, "meta.political_participation_biv.csv"))
participation.effects     <- read_if(paste0(base_path_effects, "5.results_political_participation_models.csv"))

# Leadership
meta.leadership           <- read_if(paste0(base_path_meta, "meta.leadership.csv"))
meta.leadership.random    <- read_if(paste0(base_path_meta, "meta.leadership_random.csv"))
meta.leadership_biv       <- read_if(paste0(base_path_meta, "meta.leadership_biv.csv"))
leadership.effects        <- read_if(paste0(base_path_effects, "2.results_leadership_models.csv"))

# Interest
meta.interest             <- read_if(paste0(base_path_meta, "meta.political_interest.csv"))
meta.interest.random      <- read_if(paste0(base_path_meta, "meta.political_interest_random.csv"))
meta.interest_biv         <- read_if(paste0(base_path_meta, "meta.political_interest_biv.csv"))
interest.effects          <- read_if(paste0(base_path_effects, "6.results_political_interest_models.csv"))

# Turnout
meta.turnout              <- read_if(paste0(base_path_meta, "meta.voting.csv"))
meta.turnout.random       <- read_if(paste0(base_path_meta, "meta.voting_random.csv"))
meta.turnout_biv          <- read_if(paste0(base_path_meta, "meta.voting_biv.csv"))
turnout.effects           <- read_if(paste0(base_path_effects, "8.results_voting_models.csv"))

# Trust
meta.trust                <- read_if(paste0(base_path_meta, "meta.generalized_trust.csv"))
meta.trust.random         <- read_if(paste0(base_path_meta, "meta.generalized_trust_random.csv"))
meta.trust_biv            <- read_if(paste0(base_path_meta, "meta.generalized_trust_biv.csv"))
trust.effects             <- read_if(paste0(base_path_effects, "3.results_generalized_trust_models.csv"))

# Altruism
meta.altruism             <- read_if(paste0(base_path_meta, "meta.altruism.csv"))
meta.altruism.random      <- read_if(paste0(base_path_meta, "meta.altruism_random.csv"))
meta.altruism_biv         <- read_if(paste0(base_path_meta, "meta.altruism_biv.csv"))
altruism.effects          <- read_if(paste0(base_path_effects, "1.results_altruism_models.csv"))

# Normative games
meta.normativegames        <- read_if(paste0(base_path_meta, "meta.normative.csv"))
meta.normativegames_biv    <- read_if(paste0(base_path_meta, "meta.normative_biv.csv"))
normativegames.effects     <- read_if(paste0(base_path_effects, "4.results_normative_models.csv"))

# Stack effects
re3_estimates <- bind_rows(
  groups.effects, participation.effects, leadership.effects, interest.effects,
  turnout.effects, trust.effects, altruism.effects, normativegames.effects
)

## ======================== Plot blocks ========================

## 1) Social group participation
orig_N <- n_distinct(meta.groups$authoryear); orig_k <- nrow(meta.groups)
qe_N   <- if (is.null(meta.groups.random)) 0 else n_distinct(meta.groups.random$authoryear)
qe_k   <- if (is.null(meta.groups.random)) 0 else nrow(meta.groups.random)

plot_data_groups <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Social group participation")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model == "Full - All" ~ 1.012,
      sample_model == "Full - Quasi-experimental" ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model == "Full - All" ~ col_full,
      sample_model == "Full - Quasi-experimental" ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = case_when(
      sample_model == "Full - Quasi-experimental" ~ sh_qe,
      TRUE ~ sh_full
    ),
    lab = lbl(coef, pval)
  ) %>% drop_na(y)

groups.plot <- ggplot(meta.groups, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.groups_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_point(data = plot_data_groups, aes(x = coef, y = y, color = col), shape = plot_data_groups$shp, size = 2, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = plot_data_groups, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .8) +
  geom_text(data = plot_data_groups, aes(x = coef, y = y + .005, label = lab, color = col), size = 3, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Social\ngroup\nparticipation") +
  theme_void() + theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
                       plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", orig_N, "\nk = ", orig_k), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qe_N,   "\nk = ", qe_k ), size = 3)

## 2) Political participation
orig_N <- n_distinct(meta.participation$authoryear); orig_k <- nrow(meta.participation)
qe_N   <- if (is.null(meta.participation.random)) 0 else n_distinct(meta.participation.random$authoryear)
qe_k   <- if (is.null(meta.participation.random)) 0 else nrow(meta.participation.random)

plot_data_participation <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Political participation")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model == "Full - All" ~ 1.012,
      sample_model == "Full - Quasi-experimental" ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model == "Full - All" ~ col_full,
      sample_model == "Full - Quasi-experimental" ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model == "Full - Quasi-experimental", sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% drop_na(y)

participation.plot <- ggplot(meta.participation, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.participation_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_point(data = plot_data_participation, aes(x = coef, y = y, color = col), shape = plot_data_participation$shp, size = 2, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = plot_data_participation, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .8) +
  geom_text(data = plot_data_participation, aes(x = coef, y = y + .005, label = lab, color = col), size = 2.5, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Political\nparticipation") +
  theme_void() + theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
                       plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", orig_N, "\nk = ", orig_k), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qe_N,   "\nk = ", qe_k ), size = 3)

## 3) Community leadership
orig_N <- n_distinct(meta.leadership$authoryear); orig_k <- nrow(meta.leadership)
qe_N   <- if (is.null(meta.leadership.random)) 0 else n_distinct(meta.leadership.random$authoryear)
qe_k   <- if (is.null(meta.leadership.random)) 0 else nrow(meta.leadership.random)

plot_data_leadership <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Leadership")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model == "Full - All" ~ 1.012,
      sample_model == "Full - Quasi-experimental" ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model == "Full - All" ~ col_full,
      sample_model == "Full - Quasi-experimental" ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model == "Full - Quasi-experimental", sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% drop_na(y)

leadership.plot <- ggplot(meta.leadership, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.leadership_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_point(data = plot_data_leadership, aes(x = coef, y = y, color = col), shape = plot_data_leadership$shp, size = 2, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = plot_data_leadership, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .8) +
  geom_text(data = plot_data_leadership, aes(x = coef, y = y + .005, label = lab, color = col), size = 2.5, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Community\nleadership") +
  theme_void() + theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
                       plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", orig_N, "\nk = ", orig_k), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qe_N,   "\nk = ", qe_k ), size = 3)

## 4) Political interest / knowledge
orig_N <- n_distinct(meta.interest$authoryear); orig_k <- nrow(meta.interest)
qe_N   <- if (is.null(meta.interest.random)) 0 else n_distinct(meta.interest.random$authoryear)
qe_k   <- if (is.null(meta.interest.random)) 0 else nrow(meta.interest.random)

plot_data_interest <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Political interest")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model == "Full - All" ~ 1.012,
      sample_model == "Full - Quasi-experimental" ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model == "Full - All" ~ col_full,
      sample_model == "Full - Quasi-experimental" ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model == "Full - Quasi-experimental", sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% drop_na(y)

interest.plot <- ggplot(meta.interest, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.interest_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_point(data = plot_data_interest, aes(x = coef, y = y, color = col), shape = plot_data_interest$shp, size = 2, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = plot_data_interest, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .8) +
  geom_text(data = plot_data_interest, aes(x = coef, y = y + .005, label = lab, color = col), size = 2.5, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Political\ninterest /\nknowledge") +
  theme_void() + theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
                       plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", orig_N, "\nk = ", orig_k), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qe_N,   "\nk = ", qe_k ), size = 3)

## 5) Voting
orig_N <- n_distinct(meta.turnout$authoryear); orig_k <- nrow(meta.turnout)
qe_N   <- if (is.null(meta.turnout.random)) 0 else n_distinct(meta.turnout.random$authoryear)
qe_k   <- if (is.null(meta.turnout.random)) 0 else nrow(meta.turnout.random)

plot_data_turnout <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Voting")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model == "Full - All" ~ 1.012,
      sample_model == "Full - Quasi-experimental" ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model == "Full - All" ~ col_full,
      sample_model == "Full - Quasi-experimental" ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model == "Full - Quasi-experimental", sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% drop_na(y)

turnout.plot <- ggplot(meta.turnout, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.turnout_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_point(data = plot_data_turnout, aes(x = coef, y = y, color = col), shape = plot_data_turnout$shp, size = 2, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = plot_data_turnout, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .8) +
  geom_text(data = plot_data_turnout, aes(x = coef, y = y + .005, label = lab, color = col), size = 2.5, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Voting") +
  theme_void() + theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
                       plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", orig_N, "\nk = ", orig_k), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qe_N,   "\nk = ", qe_k ), size = 3)

## 6) Generalized trust
orig_N <- n_distinct(meta.trust$authoryear); orig_k <- nrow(meta.trust)
qe_N   <- if (is.null(meta.trust.random)) 0 else n_distinct(meta.trust.random$authoryear)
qe_k   <- if (is.null(meta.trust.random)) 0 else nrow(meta.trust.random)

plot_data_trust <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Generalized Trust")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model == "Full - All" ~ 1.012,
      sample_model == "Full - Quasi-experimental" ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model == "Full - All" ~ col_full,
      sample_model == "Full - Quasi-experimental" ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model == "Full - Quasi-experimental", sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% drop_na(y)

trust.plot <- ggplot(meta.trust, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.trust_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_point(data = plot_data_trust, aes(x = coef, y = y, color = col), shape = plot_data_trust$shp, size = 2, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = plot_data_trust, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .8) +
  geom_text(data = plot_data_trust, aes(x = coef, y = y + .005, label = lab, color = col), size = 2.5, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Generalized\ntrust") +
  theme_void() + theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
                       plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", orig_N, "\nk = ", orig_k), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qe_N,   "\nk = ", qe_k ), size = 3)

## 7) Altruism
orig_N <- n_distinct(meta.altruism$authoryear); orig_k <- nrow(meta.altruism)
qe_N   <- if (is.null(meta.altruism.random)) 0 else n_distinct(meta.altruism.random$authoryear)
qe_k   <- if (is.null(meta.altruism.random)) 0 else nrow(meta.altruism.random)

plot_data_altruism <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Altruism")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model == "Full - All" ~ 1.012,
      sample_model == "Full - Quasi-experimental" ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model == "Full - All" ~ col_full,
      sample_model == "Full - Quasi-experimental" ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model == "Full - Quasi-experimental", sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% drop_na(y)

altruism.plot <- ggplot(meta.altruism, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.altruism_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_point(data = plot_data_altruism, aes(x = coef, y = y, color = col), shape = plot_data_altruism$shp, size = 2, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = plot_data_altruism, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .8) +
  geom_text(data = plot_data_altruism, aes(x = coef, y = y + .005, label = lab, color = col), size = 2.5, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Altruism") +
  theme_void() + theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
                       plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", orig_N, "\nk = ", orig_k), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qe_N,   "\nk = ", qe_k ), size = 3)

## 8) Normative prosociality
orig_N <- n_distinct(meta.normativegames$authoryear); orig_k <- nrow(meta.normativegames)
qe_N   <- 0
qe_k   <- 0

plot_data_normative <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Normative Prosociality")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model == "Full - All" ~ 1.012,
      sample_model == "Full - Quasi-experimental" ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model == "Full - All" ~ col_full,
      sample_model == "Full - Quasi-experimental" ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model == "Full - Quasi-experimental", sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% drop_na(y)

normative.plot <- ggplot(meta.normativegames, aes(x = coef, y = 1, size = 1 / se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = 0.3, width = 0, height = 0.012) +
  geom_jitter(data = meta.normativegames_biv, aes(x = coef, y = 1, size = 1 / se),
              shape = 4, stroke = 1, color = "gray70", alpha = 0.2, width = 0, height = 0.012) +
  geom_point(data = plot_data_normative, aes(x = coef, y = y, color = col), shape = plot_data_normative$shp, size = 2, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = plot_data_normative, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = 0.8) +
  geom_text(data = plot_data_normative, aes(x = coef, y = y + 0.005, label = lab, color = col), size = 2.5, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "Effect size", y = "Normative\nprosociality") +
  theme_minimal() +
  theme(
    axis.title.x = element_text(size = 10, face = "bold"),
    axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
    axis.text.y = element_blank(), axis.ticks.y = element_blank(),
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    plot.margin = margin(0, 0, 0, 0), legend.position = "none"
  ) +
  xlim(c(-0.4, 0.4)) + ylim(c(0.985, 1.025)) +
  ggplot2::annotate("text", x = 0.30, y = 1.000, label = paste0("All:\nN = ", orig_N, "\nk = ", orig_k), size = 3) +
  ggplot2::annotate("text", x = 0.40, y = 1.000, label = paste0("QE:\nN = ", qe_N,   "\nk = ", qe_k ), size = 3)

## ------------------------ Legend and export ------------------------
legend_plot <- ggplot() +
  xlim(0, 1) + ylim(0, 1) + theme_void() +
  ggplot2::annotate("point", x = 0.70, y = 0.90, shape = sh_full, size = 3, color = col_full, fill = col_full) +
  ggplot2::annotate("text",  x = 0.73, y = 0.90, label = "Full models", hjust = 0, vjust = 0.5, size = 3) +
  ggplot2::annotate("point", x = 0.70, y = 0.70, shape = sh_biv,  size = 3, color = col_biv,  fill = col_biv) +
  ggplot2::annotate("text",  x = 0.73, y = 0.70, label = "Bivariate models", hjust = 0, vjust = 0.5, size = 3) +
  ggplot2::annotate("point", x = 0.70, y = 0.50, shape = sh_qe,   size = 3, color = col_qe,  fill = col_qe) +
  ggplot2::annotate("text",  x = 0.73, y = 0.50, label = "Quasi experimental models", hjust = 0, vjust = 0.5, size = 3) +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  labs(x = "", y = "")

combined_plot <- (
    groups.plot /
    participation.plot /
    altruism.plot /
    leadership.plot /
    interest.plot /
    turnout.plot /
    trust.plot /
    normative.plot /
    legend_plot
) & theme(axis.title.x = element_text(face = "bold")) & plot_layout(guides = "collect")

#out_path <- "/Users/jbs548/Library/CloudStorage/Dropbox/Joan Barcelo/Present/NYUAD Assistant Professor/Research/Papers/Work in Progress/Meta_analysis/Analysis/Replication_materials_APSR/Figures/"
#pdf(paste0(out_path, "Figure2.pdf"),
#    height = 13, width = 8.4, paper = "letter")
#print(combined_plot)
#grid.draw(linesGrob(x = unit(c(0.528, 0.528), "npc"),
#                    y = unit(c(0.16, 0.97), "npc"),
#                    gp = gpar(col = "darkgray", alpha = 0.4)))
#dev.off()